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The transition from “few to many” lias recently been probed experimentally in an ultracold 
harmonically confined one-dimensional lithium gas, in which a single impurity atom interacts with 
a background gas consisting of one, two, or more identical fermions [A. N. Wenz et al., Science 
342 , 457 (2013)]. For repulsive interactions between the background or majority atoms and the 
impurity, the interaction energy for relatively moderate system sizes was analyzed and found to 
converge toward the corresponding expression for an infinitely large Fermi gas. Motivated by these 
experimental results, we apply perturbative techniques to determine the interaction energy for weak 
and strong coupling strengths and derive approximate descriptions for the interaction energy for 
repulsive interactions with varying strength between the impurity and the majority atoms and any 
number of majority atoms. 

PACS numbers: 


I. INTRODUCTION 

One-dimensional Bose and Fermi systems with contact 
interactions have been studied for many decades now, 
especially in the regime where the systems obey peri¬ 
odic boundary conditions EH- A large fraction of the 
eigenstates can be thought of as corresponding to gas- 
like states. A second subset of eigenstates corresponds 
to self-bound droplet-like states. These states maintain 
their bound state character in the absence of periodic 
boundary conditions, i.e., in free space. In many cases, 
both the gas-like and droplet-like states can be obtained 
analytically via the Bethe ansatz. The Bethe ansatz takes 
advantage of the fact that the zero-range nature of the 
interactions, combined with the fact that particles in one 
dimension have to pass through each other to exchange 
positions, allows one to identify constants of motion. The 
solutions can then be derived in terms of these constants 
of motion. A closely related aspect is that a variety of 
one-dimensional systems with two-body contact interac¬ 
tions are integrable jl], @] . 

The solution of the homogeneous system can be ap¬ 
plied to one-dimensional systems under spatially vary¬ 
ing external confinement via the local density approxi¬ 
mation [M3- This approximation typically provides a 
highly accurate description for a large number of par¬ 
ticles but not necessarily for a small number of parti¬ 
cles. It is thus desirable to derive more accurate descrip¬ 
tions for small one-dimensional systems with two-body 
delta-function interactions under external confinement. 
Unfortunately, extensions of the Bethe ansatz to inho¬ 
mogeneous systems are, in general, not known. This 
can be understood intuitively by realizing that the rel¬ 
ative two-body momentum in inhomogeneous systems is 
not conserved due to the presence of the spatially vary¬ 
ing confinement. Correspondingly, harmonically trapped 
one-dimensional few-body systems have been treated nu¬ 
merically by various techniques EUHHI3- 

In this work, we apply standard Raleigh-Schrodinger 


perturbation theory to harmonically confined systems 
and derive approximate solutions whose accuracy can 
be improved systematically by considering successively 
higher orders in the expansion in the small parame¬ 
ter. We focus on one-dimensional Fermi gases with a 
single impurity under external harmonic confinement. 
This system is of particular interest since it has been 
realized experimentally in Jochim’s cold atom labora¬ 
tory EM E3|- I n the experiments, the impurity is a 
lithium atom that occupies a hyperfine state different 
from the hyperfine state that the majority atoms oc¬ 
cupy. The trapping geometry is highly-elongated and ef¬ 
fectively one-dimensional. We will show that our pertur¬ 
bative results enable us to calculate the energy of the up¬ 
per branch, which has been studied experimentally, with 
fairly good accuracy for all N over a wide range of cou¬ 
pling strengths. In addition, our results provide bounds 
on the energies in the weakly- and strongly-interacting 
regimes. These bounds can, e.g., be used to assess the 
accuracy of numerical solutions. 

The remainder of this paper is organized as follows. 
Section m introduces the system Hamiltonian and nota¬ 
tion. Section m summarizes our perturbative results. 
The perturbative results are analyzed in Secs. HVIandlVl 
Finally, Sec. IVII concludes. 


II. SYSTEM HAMILTONIAN 

We consider a single impurity immersed in a one- 
dimensional Fermi gas that consists of N identical mass 
to fermions. The mass of the impurity is equal to that 
of the majority or background particles. The impurity, 
with position coordinate zq, interacts with the majority 
particles, with position coordinates Zj (j = l,--- , IV), 
through a zero-range two-body potential with strength 
9, 


Vzbizjo) = gS(zj - Zq), 


(1) 
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where zjq = Zj —Zo. The Hamiltonian H for the harmon¬ 
ically confined (N, 1) system then reads 

N N 

H = J2 tfhofe) + H ho {zo) + V 2h ( Zj o), (2) 

i=l 2=1 

where the single particle harmonic oscillator Hamiltonian 
#ho(~) is given by 

H ^ )= -L^ + i mu2z2 ’ <3) 

here, ui denotes the angular trapping frequency. The 
delta-function interactions in Eq. 0 can be replaced by 
a set of boundary conditions on the many-body wave 
function ’f'(zo) %i, ■■ ■ , 2jv), 


5 ^ 


dz 


jO 


Zj O ->0+ 


dzjo 


Zjo —>0 




( 4 ) 


where the limits Zjo —»• 0 + , Zjo —t 0 , and Zj o —> 0 
are taken while keeping the other N coordinates, i.e., 
zi, ■ ■ ■ , Zj-!,z j+ i, ■■■ ,z N and (zj + z 0 )/2, fixed. 

In the following, we determine the eigenenergies E{N ) 
of the Hamiltonian H for various N. Throughout, we 
restrict ourselves to the so-called upper branch. This 
branch can be populated by preparing the system in the 
non-interacting limit (g —> 0 + ) and by then adiabati- 
cally first increasing g to large positive values, then con¬ 
tinuing across the confinement-induced resonance [20j to 
infinitely negative g values and finally increasing g to 
small negative values. Solid, dotted, and dashed lines 
in Fig. [Ha) show the energy of the upper branch for 
N = 1 [2lJ, 2, and 3 fl-H. I7H.!22| . respectively, as a func¬ 
tion of — 1 /g. For all IV, the energy increases monotoni- 
cally as a function of increasing — 1 /g. The upper branch 
corresponds to the ground state of the model Hamilto¬ 
nian when g is positive but not when g is negative. For 
negative g, the model Hamiltonian supports molecular- 
like bound states. In real cold atom systems, energet¬ 
ically lower-lying molecular states exist even for posi¬ 
tive g. However, it has been demonstrated experimen¬ 
tally M that the upper branch can be populated with 
reasonably high fidelity for positive .<?, motivating us—as 
well as others am fl3l Ha . ll7L — to investigate 

the properties of the upper branch within a stationary 
zero-temperature quantum mechanics framework. Since 
decay to states with molecular character can lead to sig¬ 
nificant depopulation of the upper branch for negative <7, 
our primary focus in the following lies on the positive g 
portion of the upper branch. 

For g = 0 + , the energy of the upper branch is equal 
to E n [(N) = ( N 2 + l)Hu}/ 2 . We write the energy E(N) 
of the upper branch in terms of the interaction energy 
<N), 


E(N) = Eni(N) + e(N). (5) 


Solid, dotted, and dashed lines in Fig. [ljb) show the 
interaction energies, normalized by the energy Ep(N), 



FIG. 1: (Color online) (a) Solid, dotted, and dashed lines 
show the energy of the upper branch for N = 1, 2, and 3 
as a function of —1 /g. The energies of the (1,1) system are 
obtained by solving the transcendental equation derived in 
Ref. [2ll] . The energies of the (2,1) and (3, 1) system are 
taken from Refs. [0 . ITH . EH ] . (b) Solid, dotted, and dashed 
lines show the interaction energy e(N), normalized by the 
Fermi energy Ef(N), for N — 1, 2, and 3, respectively, as 
a function of —1 /g. The harmonic oscillator length dho is 
defined in Eq. USD- 


for systems with iV = 1, 2, and 3 majority particles. The 
energy Ep(N) is directly proportional to the number of 
majority particles, 


Ep(N) = NHuj. (6) 

Figure [U]b) shows that the normalized interaction en¬ 
ergy depends relatively weakly on the number of parti¬ 
cles. Independent of N, we have e(N) = 0 for g = 0 + and 
e(N) = Ep(N) for |g| = 00 . As can be read off Figs. [U a) 
andQJb), the energy increase of the upper branch is the 
same on the positive g side as it is on the negative g side, 
indicating that e(N) approaches 2Ep(N) in the g = 0“ 
limit for N = 1 — 3. We refer to Ep(N) as the Fermi 
energy of the majority particles. It should be noted, 
however, that the “exact” Fermi energy of the majority 
particles is Ef{N) — Huj/2, i.e., Ep(N) corresponds to the 
leading order term of the Fermi energy of the majority 
particles in the large N limit. 

One of the main goals of this paper is to derive ex¬ 
pansions for the interaction energy of the upper branch 
around g = 0 + and \g\ = 00 using standard Raleigh- 
Schrodinger perturbation theory for any N, i.e., for 
N = 1, ■ ■ • , 00 . To this end, we express the interaction 
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energies e^ 0+ - ) and e(l°°l) in the vicinity of g = 0 + and 
|g| = oo, respectively, in a power series of the dimension¬ 
less interaction parameter 7 (for \g\ small) or in a power 
series of 7 -1 (for |g| large) [271 ] . 


] (N) = 


B^ k \N) 7 * 


Lfc=i 


E F (N) + 0(7 k 


K + l 


) (7) 


and 


1+ J2 C {k \N)^~ 


k =1 


e (l°°l)(7V) = 

E F (N)+0( 7 - (fem -+ 1) ), ( 8 ) 


where the dimensionless interaction parameter 7 is given 

by 


with 


7 = 


7T 9 
V2N fuoa ho ’ 



(9) 


( 10 ) 


denoting the harmonic oscillator length. As we will show 
below, the scaling of the interaction energy by E F (N) 
ensures a smooth connection between the energy shifts 
for finite and infinite N. In Eqs. 0 -®, the dimension¬ 
less /cth-order perturbation theory coefficients B( k \N) 
and ( N ) depend on N and will be determined in the 
next section. 


III. PERTURBATIVE RESULTS 


TABLE I: Coefficients B^ k \N) for various (IV, fc) combina¬ 
tions. The numbers in parenthesis denote the uncertainty 
that arises from evaluating the perturbation theory sums with 
a finite energy cutoff. The numbers without errorbars have 


been rounded. 

k = 1 k = 2 k = 3 

N= 1 

0.179587 -0.0223551 0.00179230 

N =2 

0.190481 -0.0239838 0.00179523 

N =3 

0.194409 -0.0244852 0.00177603(1) 

N=4 

0.196423 -0.0247210 0.0017627(1) 

N =5 

0.197647 -0.0248563 0.0017535(1) 

N =6 

0.198469 -0.0249435 0.0017470(1) 

N=7 

0.199059 -0.0250042 

N =8 

0.199503 -0.0250488 

N=9 

0.199849 -0.0250828 

N =10 

0.200126 -0.0251096 

N=n 

0.200353 -0.0251313 

N =12 

0.200543 -0.0251491 

N =00 

0.202642 -0.0253303 0.00171100 


TABLE II: Coefficients C^ k \N) for various ( N,k ) combina¬ 
tions. The numbers in parenthesis denote the uncertainty that 
arises from evaluating the perturbation theory sums with a fi¬ 
nite energy cutoff. The numbers without errorbars have been 
rounded. 



k = 1 k — 2 k = 3 

N= 1 

-3.54491 3.85603 34.3007 

N =2 

-3.17245 2.41904(1) 25.38(2) 

N =3 

-3.02854 1.8142(2) 23.78(8) 

N =4 

-2.95040 

N =5 

-2.90081 

N =6 

-2.86634 

N =7 

-2.84091 

N= 8 

-2.82133 

N =9 

-2.80578 

N =00 

-2.66667 0 21.0552 


N —^ 00 limit: The impurity problem for the homo¬ 
geneous system with positive 7 was solved by McGuire 
in 1965 [23]. Within the local density approximation the 
Fermi wave vector is replaced by the wave vector at the 
trap center such that the interaction energy of the ground 
state for the harmonically trapped system with N —> 00 
becomes 0 


e(oo) 

E f { 00 ) 


i-l + f-X + HTUctoi 

4 \27t 7 J V 27t > 


(£) 


.(ii) 


Expanding Eq. CEB around 7 = 0 + and | 7 | = 00 , respec¬ 
tively, B^ k \ 00 ) and C^ k \ 00 ) can be obtained for k = 
1,2, • ■ ■. We find B^( 00 ) = 2/ 7 r 2 , B^ ( 00 ) = — 1 /( 47 t 2 ), 
B< 3 )( 00 ) = 1/(6tt 4 ), C'W(oo) = -8/3, C< 2 >( 00 ) = 0, 
and C'( 3 )(oo) = 327 t 2 /15. The numerical values of these 
coefficients are summarized in Tables mu It is readily 
shown that the small and large 7 series, Eqs. 0 and 
®, converge for 7 < 27r and 7 -1 < (27r) _1 , respectively. 
Table M shows that C^ 2 ^(oo) vanishes. We will return 
to this finding when we discuss the ^-dependence of the 
C( k \N) coefficients. 


(1,1) system: The eigenenergies of the harmonically 
trapped ( 1 , 1 ) system can be obtained for any 7 by solv¬ 
ing a simple transcendental equation f2lj . Expanding 
the transcendental equation around the known eigenen¬ 
ergies for small and large 7 , one obtains power series in 
the interaction energy. Inverting these series, one obtains 
analytical expressions for the B^ k \ 1 ) and C^ k \ 1 ) coeffi¬ 
cients. We find B^\ 1) = 7 r -3 / 2 , B < - 2 \ 1) = — ln( 2 )/ 7 r 3 , 
£(3)(1) = -[tt 2 - 9ln(4) 2 ]/(247r 9 / 2 ), CW(1) = -27T 1 / 2 , 
C( 2 ) (1) = - 47 r[ln( 2 ) - 1], and C< 3 )(1) = 7 r 3 / 2 [ 7 r 2 - 24 - 
9(ln(4) — 4) ln(4)]/3. The numerical values of these coef¬ 
ficients are summarized in Tables H1ITI1 As in the N —> 00 
case, the small and large 7 series for N = 1, Eqs. 0 and 
®, have a finite radius of convergence. Employing the 
techniques of Ref. [29[ , we find—using up to 50 expansion 
coefficients—that the small and large 7 series converge 
for | 7 | < 1.0745(2) x 2 tt and H ” 1 < [1.0745(2) x 27 t]-\ 
respectively. Our result for the convergence of the small 
7 series is consistent with what is reported in the litera¬ 
ture [30| . 

Weakly-repulsive (N, 1) system, N = 2,3, • • • : To treat 
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the weakly-interacting system with finite TV, TV > 1, we 
rewrite the system Hamiltonian in second quantization 
and expand the field operators for the majority particles 
and the impurity in terms of single particle harmonic 
oscillator states (see, e.g., Ref. [Ml). The interaction 
matrix elements can be evaluated analytically and the 
first-order perturbation theory treatment for positive g 
yields 


2y / TVT(l/2 + TV) 

7 r 2 TV! ' ( ' 

The first-order energy shift may be interpreted 
as the leading-order mean-field shift. We find 
lirn/v^oo (TV) = 2 / 7 r 2 , which agrees with the coef¬ 
ficient obtained by expanding Eq. m- The evalua¬ 
tion of the second-order energy shift involves the eval¬ 
uation of infinite sums. We find, as expected, that these 
sums converge. The reason is that the one-dimensional 
delta-function interaction does not, unlike two- or three- 
dimensional delta-function interactions [32] , [33| , require 
any regularization if used in standard perturbation the¬ 
ory approaches. We did not find a compact analytical ex¬ 
pression applicable to all TV for the second-order energy 
shift. For TV = 1 and 2, we have B^ 2 \ 1) = — ln(2)/7r 3 
and f?( 2 )(2) = [-9 + 6 v / 3 + 31n(2 + V3) - 12 ln(2)]/(47r 3 ). 
For larger TV, the expressions are lengthy. The numer¬ 
ical values for TV < 12 are listed in Table HI Table U 
also summarizes the numerically determined values for 
the third-order coefficients (N) for TV = 2 — 6. The 
ij( 3 ) (TV) coefficient increases slightly as TV changes from 
1 to 2, and then decreases monotonically as TV increases 
further. The numerically determined B^ (TV) coefficients 
for TV = 2 — 6 approach the TV = 00 coefficient smoothly 
if plotted as a function of 1/TV. 

Strongly-interacting (TV, 1) system, TV = 2, 3, • ■ •: The 
strongly-interacting regime has been treated perturba- 
tively at leading order, i.e., at order 1 / 7 , for TV < 8 [2(| 
(note, though, that only the coefficients for TV < 4 were 
reported explicitly, i.e., in equation or numerical form). 
To derive these results, the two-body interaction for large 
| < 7 1 is modeled by imposing the two-body boundary condi¬ 
tion on the many-body wave function when the distance 
between the unlike particles approaches zero [23| [26|, [34| . 
Since the ground state eigenenergy for |y| = 00 is degen¬ 
erate, the perturbation shift is obtained by diagonalizing 
the perturbation matrix constructed using the degener¬ 
ate states for g = 00 . For the many-body states 4/ Q and 
T. 3 , the perturbation matrix element V a p reads [23| 



^a/3 = - 


h A 

m 2 g 


N 

E 

2=1 • 


3** 


dzj 1 


•j 0 






dz 


•j 0 


z io -s-0+ ® Z jO 

<9T, 


2jo-s-0+ 


S(z j0 ) X 

Zj 0 —>-0“ J 

dzgd^i • • • dzjv(13) 


ZjQ —>-0 


These matrix elements are closely related to the bound¬ 
ary condition representation of the one-dimensional odd- 
parity pseudo-potential [H, [36|. We evaluate the inte¬ 
grals in Eq. m analytically for TV = 1 — 4. The analyt¬ 
ical results for TV = 1 and 2 read (7^(1) = —2^fn and 
C7)(2) = —y / 7r/2(81/32). The analytical expressions 
for TV = 3 and 4 are lengthy and not reported here [37j . 
For larger TV, we perform all but one integration for each 
of the perturbation matrix elements analytically. The re¬ 
sulting numerically determined energy shifts are accurate 
to more than 10 digits. Table HT1 summarizes the numeri¬ 
cal values for the coefficient C W (TV) for TV < 9 obtained 
by us. The extension to larger TV is, although tedious, 
possible in principle. 

To determine the energy shift proportional to 7 -2 , 
we use second-order perturbation theory. Reference [MJ 
pointed out that the second-order perturbation theory 
energy shift of the (1,1) system diverges, thus requiring 
regularization. Analogous divergencies arise in the per¬ 
turbative treatment of one-dimensional single-component 
Fermi gases with generalized delta-function interactions 
(see, e.g., Ref. [35|) and that of one-dimensional Bose 
gases with effective range dependent zero-range interac¬ 
tions. In the following, we discuss the impurity problem 
with TV = 2 and 3. To evaluate the second-order energy 
shifts, we need to know the complete set of eigenstates 
of the (2,1) and (3,1) systems with |g| = 00 . For the 

(2.1) system, we use the analytical wave functions from 
Ref. [39[ and evaluate the integrals analytically. For the 

(3.1) system, we derive compact forms for the eigenstates 
using spherical coordinates and evaluate the relevant in¬ 
tegrals analytically. We then evaluate the second-order 
perturbation theory sums numerically, imposing an en¬ 
ergy cutoff on the relative energy of the intermediate (or 
virtual) states that are being summed over. The second- 
order energy shift is found to contain powerlaw divergen¬ 
cies in the energy cutoff. These divergencies are canceled 
through the introduction of a counterterm and the con¬ 
stant (and physically meaningful) part is extracted with 
high precision by a regularization scheme similar to that 
developed for harmonically trapped bosons [4(|. Table HD 
reports the resulting second-order perturbation theory 
coefficients with errorbars. Our perturbative coefficients 
are consistent with the coefficients obtained by fitting 
the (2,1) and (3,1) energies reported in Refs, [ljl, |22| to 
a polynomial in 7 _1 . For the TV = 2 and 3 systems, we 
extend the above treatment to the third order (see Ta- 
ble[TT|. These third-order calculations require the evalua¬ 
tion of matrix elements V a p between excited states. Since 
the third-order perturbation expression is more involved 
than the second-order perturbation expression, our third- 
order result has a larger errorbar than our second-order 
result [4l| . 

The calculations of the second- and third-order energy 
shifts can, in principle, be extended to larger TV. To do 
so, two challenges need to be overcome. First, an efficient 
method to generate the complete set of eigenstates at 
| g | = 00 has to be devised. Second, an efficient scheme for 
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TABLE III: Fitting coefficients b^ for k = 2 and 3. For k = 2 
and 3, we used j max = 6 and 4, respectively. 



k = 2 

k = 3 

o 

II 

■c-s 

-0.0253304 

1.71100 x 10 -3 

j=l 

0.0019591 

2.23905 x 10 -4 

3=2 

0.0033477 

6.59881 x 10" 6 

J=3 

-0.0116972 

-3.67025 x 10" 4 

<-o. <-o. 

II II II 

Q Ol Ji. 

0.0379414 

-0.0684440 

0.0486346 

2.64073 x 10 -4 


evaluating the matrix elements and infinite perturbation 
theory sums has to be developed. This is not pursued 
here. 

IV. FITTING THE B (k) {N) AND C W {N) 
COEFFICIENTS 

Tables HHTT] suggest that the coefficients B^ ( TV ) and 
C^ k ) (TV) change, for fixed k , smoothly with TV. This mo¬ 
tivates us to write 

imax / i \ j 

= (14) 

j=0 ' ' 

and 

Jmax / -I \ J 

C'W(TV) = £cf ) (-) . (15) 

i=o ^ ' 

It should be noted that the expressions ©D-GSD reduce 
to bg k ' 1 and c g k \ respectively, in the TV — > oo limit. In the 
following, the parameters b^ and are obtained by 
fitting the coefficients B^ k \N) and C^(N) for fixed k. 

We start with B ( ' 2 \N). We fit Eq. (fT-fl) to the B^ 2 \N) 
values for TV = 1 — 80 (the values for TV = 1 — 12 are 
reported in Table HT1). varying j max from 2 — 20. We find 
that the most reliable fit is obtained for j max = 12 — 
13. In this case, the fitting parameter b y 0 ’ differs from 
— 1/(47t 2 ) [the result obtained by expanding Eq. (fill) ] by 
less than 10 -8 . This suggests that not only the k = 
1 coefficient (see discussion above) but also the k = 2 
coefficient connects smoothly with the infinite TV result. 
Table [TTT] reports the results of our fit to the B^ (TV) 
coefficients with TV = 1 — 80 and oo by a polynomial with 

Jmax — 0. 

As mentioned earlier, the i?( 3 )(l) coefficient is slightly 
smaller than the B^( 2) coefficient. The B^(N) coeffi¬ 
cients for TV > 2, however, decrease monotonically. This 
motivates us to fit the B^ (TV) coefficients with TV = 2—6 
and oo by a polynomial with j max = 4. The fit coefficients 
are reported in Table uni It can be seen that the coef¬ 
ficient 5 q 3 ' ) agrees with the coefficient B^ 3 \ oo) reported 
in Table HI We believe that our fit provides an accurate 
description of the 6 < TV < oo coefficients. 



1/N 


FIG. 2: (Color online) Symbols show the coefficients (a) 
C (1) (TV), (b) C (2) (N), and (c) C (3) (V) as a function of 1/N. 
The solid lines show our fits with j max = 6, 3, and 3, respec¬ 
tively. 


TABLE IV: Fitting coefficients c^ k> for k — 1,2, and 3. For 
k — 1,2, and 3, we used jmax = 6, 3, and 3, respectively. 



k = 1 

k = 2 

k = 3 

3=0 

-2.66667 

0.00000 

21.05520 

3 =1 

-1.40749 

7.06739 

8.80915 

3=2 

1.78704 

-5.70589 

-5.07455 

3 =3 

-4.21746 

2.49453 

9.51090 

j=4 

7.33094 



j =5 

-7.13604 



j =6 

2.76476 




Symbols in Figs. [2|a)-[2](c) show the C^ k \N) coeffi¬ 
cients with k = 1,2, and 3, respectively, as a function of 
1/TV. Our fits to these data (see Table EU using poly¬ 
nomials with j max = 6 , 3, and 3 are shown by solid lines 
(see Table ITVl for the coefficients). 

The discussion so far has focused on the coefficients 
B^ k \N) and C^ k \N) with k = 1-3. It is, in general, not 
feasible to extend the perturbative calculations to higher 
k for arbitrary TV. However, for TV = 1 and oo, the co¬ 
efficients with larger k can be obtained readily. We find 
that |I?( fc )(l)| decreases monotonically with increasing k 
(we checked this for k < 50). The |C^ fc ^(l)| coefficient 
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increases monotonically with increasing k for k < 37; for 
k > 37, we observe small non-monotonic oscillations. For 
N = oo, we find that the B^ (oo) with k even and k > 4 
vanish while the |i?l fe )(oo)| with k odd decrease mono¬ 
tonically with increasing k (again, we checked this for 
k < 50). Similarly, the C^ k \ oo) with k even and k > 2 
vanish while the |C^ fc ^(oo)| with k odd increase monoton¬ 
ically with increasing k. Assuming a linear change with 
1/N, interpolating between B^ k \ 1) and B^ k \ oo) and be¬ 
tween C^ k \ 1) and C^ k \ oo) for k > 3 yields estimates for 
the finite N, N > 1, coefficients. While rough, these es¬ 
timates might provide a reasonable means to connect the 
weak and strong perturbation theory limits for quantities 
such as those shown in Figs. [3] and |U 

We cannot accurately estimate the radius of conver¬ 
gence of the small and large 7 expansions for 1 < N < 00 . 
However, the fact that the radius of convergence is given 
by |b| < 1.0745(2) x 2 tt for N = 1 and 7 < 27t for N = 00 
for the small 7 series and by | 7| -1 < [1.0745(2) x 27 t] - 1 
for N = 1 and 7 -1 < for N = 00 for the large 

7 series suggests two speculations: First, a convergent 
series can be found for any 7 and N. Second, the radius 
of convergence of the small 7 series is approximately 2-7 t 
for all N. Figures [3] and [4j which are discussed in the 
next section, are consistent with these speculations. 


Ci 
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V. DISCUSSION 


This section compares the perturbative energy expres¬ 
sions with the numerically determined energies of the up¬ 
per branch. Figure[]Jb) shows that the scaled interaction 
energy e(N)/E F (N) depends weakly on N if plotted as 
a function of — fojja Fo /g. The dependence on N is even 
weaker when the interaction strength is parameterized by 
7 as opposed to g. To benchmark the applicability of the 
perturbative expressions we analyze the interaction en¬ 
ergy of the system with N majority atoms by comparing 
with that of the (1,1) system. Specifically, we consider 
the quantities p(N), 


p(N) 


e(N)/E F (N) 
£ (1)/Af(1) ’ 


(16) 


and S(N, N'), 


MN N ,, = e(N)/E F (N)-e(l)/E F (l) 

1 J e(N')/E F (N')-e(l)/E F (l)- 


(17) 


For finite N, p(N) reduces to e(7V)/[./Ve(l)], i.e., p(N) 
tells one the interaction energy per particle, normalized 
by the interaction energy of the (1,1) system. The quan¬ 
tity 6(N,N') can alternatively be written as [p{N) — 
i\MN')-i}. 

Expanding Eq. m in the weakly-interacting (small 7 ) 
regime, we find 


p(N) = pW(N) + pW(N) 7 + p { ™\N) 7 2 + 0(7 3 ),(18) 


FIG. 3: (Color online) Solid lines show the quantity p(N) 
as a function of — 1/7 for (a) N = 2, (b) N = 3, and (c) 
N = 00 , respectively. For comparison, dotted, dashed, and 
dash-dotted lines show the perturbative results for p(N) ac¬ 
counting for terms up to order 7 0 , 7 1 and 7 2 , respectively, 
in the weakly-interacting regime and accounting for terms 
up to order 7 -1 , 7 _2 and 7 -3 , respectively, in the strongly- 
interacting regime. 


where the coefficients p^\N) are determined by the 
Z?(b (TV) and B^ l \ 1) with l <k + 1. Expanding Eq. (TTCJll 
in the strongly-interacting (large | 7 |) regime, we find 

P (N) = 1 + p ( s )(iV) 7- 1 +p( s) (7V ) 7 - 2 + 

p( s) (JV ) 7 - 3 + 0(7 - 4 ) 1 (19) 

where the coefficients p^(N) are determined by the 
C w (A0 and C w ( 1 ) with l < k. Solid lines in Figs.^a)- 
m c) show the quantity p{N) for N = 2, 3, and 00 , re¬ 
spectively. The solid lines are obtained using the nu¬ 
merical (2,1) and (3,1) energies and the semi-analytical 
(1,1) and (oo,l) energies. For 7 —> 0 + (i.e., for 
— 1/7 —> — 00 ), the quantity p(N) approaches the con¬ 
stant p< w) (N) = B^(N)/B^( 1) [see the horizontal dot¬ 
ted lines in Figs. El[a)-[2Kc)] , which increases monotoni¬ 
cally from 1.0607 to 1.1284 as N goes from 2 to 00 . This 
portion of the interaction energy can be interpreted as the 
mean-field contribution. Inclusion of the next order cor¬ 
rection [the p ^ (N)'y term] and the next two corrections 
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[the p l \ x \N)"f and p^\N)j 2 terms] yields the dashed 
and dash-dotted lines in Figs.[3][a)-[3|c). The dash-dotted 
lines provide a fairly accurate description of the quantity 
p(N) for —l/y < —0.4. For |y| —» 00 , the leading-order 
7 -dependent term [see the (non-horizontal) dotted lines 
in Figs.[3la)-[3](c)] increases monotonically from 0.3725 to 
0.8782 as TV changes from 2 to 00 . Inclusion of the next- 
order correction and the next two corrections yields the 
dashed and dash-dotted lines in Figs. [3][a)-[3][c). It can 
be seen that the dash-dotted lines provide a fairly accu¬ 
rate description of the quantity p(N) for — 1/7 > —0.15. 
This value is close to the expected radius of convergence 
of the interaction energy [recall, the radius of convergence 
is I /7 = (1.0745 x 27 t) - 1 ss 0.148 for the (1,1) system]. 
Combining the perturbative descriptions for small and 
large | 7 |, the expansions provide a fairly accurate de¬ 
scription of the interaction energy for the system with 
TV majority particles, normalized by that for the (1,1) 
system, over a wide range of interaction strengths 7 . 

Expanding Eq. m in the weakly-interacting regime, 
we find 

S(N, TV') = ^ w) (TV, TV') + <5< w) (TV, TV ') 7 + 

4 w) (TV,TV')7 2 + 0(7 3 ), (20) 

where the coefficients <5£ w ' ) (TV, TV') are determined by the 
B^(N), B^(N'), and B«( 1) with l < k+ 1. Expanding 
Eq. inn in the strongly-interacting regime, we find 

<5 (TV, TV') = S ( 0 s) (TV, TV') + < 5 ^ s) (TV, TV ') 7 -1 + 

4 s) (TV,TV')7- 2 +0( 7 - 3 ), (21) 

where the coefficients 8 ^ (TV, TV') are determined by the 
C'W(TV), C'W(IV'), and C^(l) with l <k+ 1. The quan¬ 
tity 8(N, TV') is shown by the solid line in Fig. 0|a) for 
(TV, TV') = (2, 00 ) and by dots in Figs. UK b) -UK c) for (3, 00 ) 
and (2,3), respectively. We observe that the quantity 
J(TV, TV') changes only slightly as — 1/7 goes from —00 
to 0; this is particularly true for S( 2,3) [see Fig. 01c)]. 
The limiting values [see the dotted lines in Figs. 0Ja)- 
Hc)] are given by 8 ( q w \n, TV') and 8^ (TV, TV'), respec¬ 
tively. Dashed lines include the next order correction in 
the weakly- and strongly-interacting regimes, and dash- 
dotted lines include the next two corrections. In the 
weakly-interacting regime, the dash-dotted lines provide 
a fairly good description of the quantity 5(TV, TV'). In the 
strongly-interacting regime, however, the validity regime 
of the perturbative expressions is quite small. For 8 ( 2 , 3), 
e.g., the expansion coefficients are 8q\2,3) = 0.7213, 
<^ s) (2,3) = 0.0694(3), and 4 S) ( 2 . 3 ) = -2.31(12), where 
the numbers in brackets denote the errorbars due to the 
uncertainties of the second- and third-order perturbation 
theory coefficients. The fact that |<$ 2 ^( 2 , 3)| » |<5^(2, 3)| 
is responsible for the turn-around of the dash-dotted line 
for large positive 7 . We note that the errorbar of the 
quantities <5(3, 00 ) and (5(2,3), obtained from the numer¬ 
ical energies [see dots in Figs. 0](b)-0|c)], is too large in 



_1_ 1 _1_1_1_u 

-1.5 -1 -0.5 0 

-l/y 


FIG. 4: (Color online) The quantity S(N, TV') as a function of 
— 1 / 7 . The solid line is for (a) TV = 2 and TV' = 00 and the 
circles are for (b) TV = 3 and N' = 00 and (c) TV = 2 and TV' = 
3. In the large 7 regime, the uncertainty of the numerically 
determined (3,1) energies leads to appreciable uncertainties 
in 8(2, 3) and 5(3, 00 ) [see the errorbars in Figs. [2b) and[2c)]. 
For comparison, dotted, dashed, and dash-dotted lines show 
the perturbative results for 5(TV, TV') accounting for terms up 
to order 7 0 , 7 1 , and y 2 , respectively, in the weakly-interacting, 
small I 7 I regime and accounting for terms up to order 7 -1 , 
7 -2 , and y~ 3 , respectively, in the strongly-interacting, large 
I 7 I regime. 


the large 7 regime to meaningfully compare with the per¬ 
turbative results. 


VI. CONCLUSION 

This paper considered the upper branch of a non¬ 
interacting harmonically trapped one-dimensional Fermi 
gas with a single impurity. Zero-range two-body contact 
interactions with strength g were assumed between the 
majority atoms and the impurity. This system consti¬ 
tutes one of the simplest mesoscopic systems accessible 
to experiment and theory. On the experimental side, it 
has been demonstrated by the Heidelberg group that the 
upper branch of the model Hamiltonian can be emulated 
reliably using ultracold atoms mm- On the theory 
















side, various numerical and_ analytical techniques have 
been applied Ifi. ITTl - fl7l [23l42(| . This paper pursued a per¬ 
turbative approach, which determined expansions of the 
energy of the upper branch in the weakly- and strongly- 
interacting regimes for various N. In the cases where 
we were not able to obtain general N expressions for a 
fixed order in the perturbative expansion, approximate 
expressions applicable to all N were obtained through 
fits. Through comparison with accurate numerical few- 
body energies, the perturbative expressions were shown 
to provide a satisfactory description for a wide range of 
interaction strengths. 

The main results of this work are: (i) We determined 
an expansion for the energy of the upper branch of a 
one-dimensional harmonically trapped Fermi gas with a 
single impurity in the weakly-repulsive regime up to or¬ 
der 7 3 , applicable to any system size, (ii) We deter¬ 
mined an expansion for the energy of the upper branch 
in the strongly-interacting regime up to order 7 , ap¬ 

plicable to any system size. While the idea to treat the 
coupling strength 1/7 as a small parameter is not new, 
our work provides an explicit demonstration that such 


a program can be carried through explicitly beyond the 
leading-order correction, (iii) The radii of convergence 
of the series were reported for N = 1 and 00 . (iv) The 
behavior of the expansion coefficients in the series in 
and 7 -fc with k > 3 was discussed, (v) The perturba¬ 
tive expressions were benchmarked and found to provide 
a reliable description over a wide range of interaction 
strengths. 

The results presented in this work can be used to cal¬ 
culate perturbative expressions for the contact and other 
observables. Moreover, the second- and third-order re¬ 
sults in the 7 -1 series allow one to assess the applicability 
regime of effective spin models [H, H 3 , |t 3 . 
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